Fluctuation-response relation in turbulent systems 
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We address the problem of measuring time-properties of Response Functions (Green functions) 
in Gaussian models (Orszag-McLaughin) and strongly non-Gaussian models (shell models for tur- 
bulence). We introduce the concept of halving time statistics to have a statistically stable tool to 
quantify the time decay of Response Functions and Generalized Response Functions of high order. 
We show numerically that in shell models for three dimensional turbulence Response Functions 
are inertial range quantities. This is a strong indication that the invariant measure describing the 
shell-velocity fluctuations is characterized by short range interactions between neighboring shells. 

O' 

O ■ I. INTRODUCTION 

<N : 

^ . Fluctuation-Response (F/R) relation plays an important role in statistical mechanics and, more generally, in systems 

H-~5 ' with chaotic dynamics. With the term F/R relation one indicates the connection between the relaxation properties of 

fT^ , a system and its response to an external perturbation. The relevance of this relation is evident: it allows to connect 
"non-equilibrium" features (i.e. response and relaxation) to "equilibrium" |l|] properties (correlation functions). As 

' -_ '\ an important example, we mention the Green-Kubo H formulas in the linear response theory which links the response 

\—A • to an external field with correlations computed at equilibrium. 

Q^ ' Consider a system whose state is given by a finite dimension vector x = {xi, ...,xn), the average hnear response 

r^ • Gl{t) = {Rl{t)) is the average response after a time t of the variable Xi to a small perturbation of the variable Xj at 

•jH ■ time t = 0. Under rather general conditions (basically one has to assume that the system is mixing) it is possible to 

f^ . show that a generalized F/R relation holds PJfI] : 



<^K*)) = (t4^)-(^*W/.W0)])' (1) 



Sx,{t) 
f— ^ , where the functions fj depend on the invariant probability distribution p(x): 

o 

i>: f fxi - ^i^^w (2) 

o. ^'^''^- dx, ■ ^^^ 
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^^ I The physical meaning of (|l]) is the following: consider a small perturbation (5x(0) — ((5a;i(0), ..., (5a;jv(0)) at t = 0; the 

r^ ' average distance {Sxi(t)) from the unperturbed values Xi{t) is 

• 1—1 ■ 

'S: (<5x.(t))=^(i?j-(i)>fe,(0). (3) 

> 

J^ For Hamiltonian systems one realizes that {^ is the usual linear response theory. If p(x) is Gaussian, one has a simple 

M relationship between the response and the correlation function: 
C^ ■ 

_ _ I 

[XiXj) — [Xi)[Xj) 

In the general case, i.e. non-Gaussian statistics, formula (|lj) gives just a qualitative information, i.e. the existence of 
a link between response and the general correlation function {xi{t)fj[:x.(0)]). 

In particular, in the most interesting cases in which p{x) is unknown, it is extremely important that F/R relation (|lj) 
exists, because it allows to control some properties of the invariant measure, p(x), in terms of the Response-Functions 
behavior. In the past, this has not always been clear, e.g. some authors claim (with qualitative arguments) that in 
fully developed turbulence there is no relation between equilibrium fluctuations and relaxation to equilibrium |12| , 
while a proper statement would limit to the non existence of the usual "Gaussian-like" F/R relation (i|). 
Response Functions have a clear phenomenological importance in many applied problems where one needs to con- 
trol and/or predict the system reaction as a function of external spatial and/or temporal perturbations. Response 
functions, also known as Green functions, play also a very important role in many theoretical attempts to attack 
non-equilibrium problems. In particular, in many analytical approach to hydrodynamical problems described by 



Navier-Stokes eqs, or models of them, Greens functions naturally show up both in perturbative schemes [BJ and 
closure-like attempts like the DIA approximation |g]. 

We stress that the F/R , in the form (|l|), is a rather general relation which does not depend too much on the 
details of the measure of the systems, e.g. both in presence or absence of an energy flux. For example, in the field 
of disordered systems the F/R had been widely studied in order to highlight non trivial relaxation aspects e.g. aging 
phenomena. In this paper we want to address the problem of F/R relation for the case of dynamical models with 
many degrees of freedom and many characteristic times. We are also interested in exploiting the F/R relation in 
models which exhibit strong departure from Gaussian statistics. We introduce a suitable new numerical method for 
measuring the characteristic times involved in the response functions. 

This new method is based on the idea to characterize the response behavior as a function of its halving time statistics 
(HTS), i.e. the time, t, necessary for the response from a typical infinitesimal perturbation to reach, say, one half of 
its initial value. 

The plan of the paper is as follows. First we investigate a dynamical system with many degrees of freedom and 
many different characteristic times where still a classical Gaussian set of F/R relations holds. The model is the so- 
called "Orszag-McLaughlin" model which is used to probe the effective improvement of halving-time-statistics with 
respect to the usual direct measurement of time decaying properties. Than, we attack the much less trivial case of 
characterizing response behavior in models for three dimensional turbulent energy cascade, i.e. shell models |g|. 
Also in the latter case, the halving time statistics will allow us to measure with good accuracy the non-trivial time 
properties of Green functions. As a result, we show that also the response function (which probes linear features of 
the dynamical evolution) is strongly affected by the non-linear inertial range physics. As a consequence, short range 
interactions between neighboring shells are thought to characterize the invariant measure describing shell-velocity 
fluctuations. 

II. NUMERICAL SIMULATIONS 

Before entering the detailed description of the results, we want to discuss a practical problem for the numerical 
computation of G*(t) = {R'j{t)). In numerical simulations, (Rjit)) is computed perturbing the variable xi at time 
t = ti with an "infinitesimal" kick of amplitude Sxj{ti\ti) — Xj{t — ti) — Xj(t = ti) = e, for e —^ 0, and the deviation 
(5x(f|ii) = x(t) — x(t) is computed integrating the two trajectories x and x up to a prescribed time t2 = ^i + ^t- At 
time t = t2 the variable Xi is again perturbed with another kick, and a new sample Sx{t) is computed and so forth. 
The procedure is repeated M <C 1 times and the mean response is then evaluated as 

r'J7+\ - ^ Y^ Sxi{tk+t\tk) , , 

In presence of chaos, the absolute value of the deviation, \Sxi(tk + t,\tk)\, typically grows exponentially with t. 
Therefore, the mean response, (RUt)), is the result of a delicate balance of terms with non-fixed sign. As a result, we 
have that the error a{t) on (RUt)) increases exponentially with t 

ait) ^ ^^ (6) 

where 7 is the generalized Lyapunov exponent of second order (greater or equal the maximum Lyapunov exponent). 
One easily understands the main problem in trying to numerically compute any response function for large times: one 
needs to control an observable which is rapidly decaying to zero with exponentially-large fluctuations. In practice, it 
turns out to be impossible to have a reliable control on the asymptotic behavior of Green functions (see next sections 
and figures therein). 

In order to avoid this trouble we propose another approach. Let us first consider only diagonal responses, i.e. response 
after a time t of the n-th variable from a perturbation of the same n-th variable at time t = 0, G^(i) = (i?"(t)). 

In this case, we claim that it is possible to have a good characterization of the main temporal properties by looking 
at the halving time statistics (HTS), that is at the probability density functions, V{t) of the time r necessary to see an 
appreciable decay of the response function: R^iT — t) — -^^"(0), with the threshold A fixed to a macroscopic value, 
say A — 1/2. In practice, one performs many response experiments, by collecting the statistics of the times necessary 
to see the response become one half of its initial value. The advantage of this HTS with respect to the more standard 
way of characterizing the mean response G," (t) with some typical time is that one docs not need to know any functional 
behavior for the averaged response and, moreover, one has also a control on the fluctuations of the characteristic times. 



i.e. the HTS integrates all times corresponding to halving events. In the following, we show that the HTS is at least 
able to reproduce with good accuracy the same results of the direct fitting procedure of the averaged response in cases 
when the classical F/R relation (Q) holds (Orszag-McLaughlin model, i.e. Gaussian statistics) and, more interesting, 
it is also able to give new hints on the F/R relations when time-intermittency and strong departure from Gaussianity 
are present (shell models). In the following wc will also discuss the cases of non-diagonal responses, (i?"j(i)), with 
n ^ m and the cases of generalized higher order responses 



A. The Orszag-McLaughlin model 

Let us consider the following model : 

— TT — ^n+lXn+2 + Xn-lXn-2 ~ -^^n+lXn-l, (') 

at 

with n = (1, 2, ..., TV), N = 20, and the periodic condition Xn+N = a;„. This model contains some of the main features 

of inviscid hydrodynamics: a) there are quadratic interactions; b) a quadratic invariant exists {E = J2n=i ^n)' '^) ^^^ 
Liouville theorem holds. For sufficiently large N the distribution of each variable a;„ is Gaussian. In this situation, 
classical F/R relationship exists for each of the n variables: self-response functions to infinitesimal perturbations are 
indistinguishable from the corresponding self-correlation functions ^J • 

We have slightly modified the system (R) in order to have variables with different characteristic times. This can be 
done, for instance, by rescaling the evolution time of each variable: 

-—— = kn{Xn+lXn+2 + Xn-lXn-2 - 2Xn+lXn-l), (8) 

where the factor fc„ is a function of the "number of identification" (e.g. site in the chain) of the variables defined as 
kn = a ■ /3", with a — 5 ■ lO^'^ and (3 = 1.7, for n ~ 1,2, ...,N/2, with the "mirror" property fc„+Ar/2 — fcAr/2+i-n- 
An immediate consequence is that the quadratic observable E is no longer invariant during the time evolution of 
the system (|^). The mean energy per mode, £'„ = {x'^) (not shown), follows a linear law with k = kn- It can be 
demonstrated that a new quadratic integral of motion exists, and this has the form: 



N 

E 



x^ 



Moreover, the a;„ variables are shown to preserve the Gaussian statistics to a good extent. Therefore, the only effect 
of the change in the original Orszag-McLaughlin system is that each variable now has its own characteristic time. 

Let us see how correlation and response functions behave for the system (g). In Figure |l| the self-correlation 
functions 

^ (^^ _ {Xn{t)Xn{0)) - {Xnf ,^„. 

and the self-response functions 

are shown. As a consequence of the preserved Gaussian statistics, F/R relation of the form (|j) holds for each of the 
variables, at least over time delays not too long. The (linear) response functions are computed as decay functions 
of single variable perturbations to infinitesimal instantaneous "kicks", averaged over a large number of simulations. 
If we conventionally define the correlation time of a variable Xn as the time delay Tc{n) after which the correlation 
function becomes lower than the value 1/2, wc find that: 

Tc{n) ^ k-^/^ (12) 



The exponent of the scahng law (|2|) can be explained with a dimensional argument, by noticing that from the mean 

energy per variable we get x^ ~ fc„, so from (g) and @ the characteristic time results to be Tc{n) ^ kn ■ We 
notice that the scaling (|l^ ) is robust with respect to the choice of the threshold value, A, i.e. it is observed even if 
the decay factor A is chosen slightly different from 1/2. 

The response time Tij(n) is defined as the time interval after which the averaged response function becomes lower 
than 1/2. We must observe that the computation of the mean response function is practically impossible after a 
certain time delay, because of exponentially growing errors. 

Last, halving times T{n) have been computed always for the variables of the system (0), using the same procedure 
as before (i.e. infinitesimal kicks). The halving time PDF's decay exponentially and, as can be shown, they all can be 
"collapsed" to the same renormalized PDF for a proper rescaling of the halving time (see below). A comparative plot 
of {T{n)), Tc{n) and Tfi{n) is shown in Figure 0. Halving times and correlation times follow the same scaling law with 
kn and, for each variable, have values very close to each other. Typical time decaying of the averaged Response, TR{n_ 
are very difficult to estimate due to high errors for slow variables. In fact, only a few points are shown in Figure 
the ones for which the mean response drops down to 1/2 fast enough, before the statistical error become too large 
(say larger than 100%). The advantage of the HTS with respect to the mean response function is that, with the same 
statistics, halving times can be computed for all variables within reasonable uncertainty, while response times. Tain), 
are generally affected by exponentially growing errors and are practically not defined when the typical relaxation time 
scale is longer than the error growth time scale. 

The numerically computed PDF's of the halving time r can be rescaled as follows: 

r^^, Pir)^{r)Pir). (13) 

We show in Figure H the overlap of some rescaled PDF's of the halving times. As a consequence, all moments of the 
T{n) PDF's have a simple scaling 
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We have found that in a Gaussian case, the response to infinitesimal perturbations can be characterized both with 
the classic mean response function and with the mean halving time technique. It is worth stressing that the HTS 
could be the only technique usable for studying relaxation to non linear perturbation in complex systems. 

B. Shell model 

Shell models for turbulent energy cascade have proved to share many statistical properties with turbulent three 
dimensional velocity fields 0,^. Let us introduces a set of wavenumber fc„ = 2"fco with n — 0,. . . ,N. The shell- 
velocity variables M„(i) must be understood as the velocity fluctuation over a distances In — k^^ . It is possible to 
write down many different sets of coupled ODEs possessing the same kinematical features necessary to mimics Navier- 
Stokes non-linear evolution. In the following we will present numerical results for a particular choice, the so-called 
Sabra model [[l6| , namely: 

[ — + vkl^ \un=i [knUn+iUn+2 + &fcn-lW„+lU* „i + (1 + h)kn-2Un-2Un-l\ + /ri, (14) 

where 6 is a free parameter, v is the molecular viscosity and /„ is an external forcing acting only at large scales, 
necessary to maintain a stationary temporal evolution. The main, strong, difference with the model discussed in the 
previous section consists in the existence of a mean energy flux from large to small scales which drives the system 
toward a strongly non-Gaussian stationary temporal evolution pq | . Shell models here discussed presents exactly the 
same qualitative difficulties of the original Navier-Stokes eqs: strong non-linearity and far from equilibrium statistical 
fluctuations. The most striking quantitative feature of the non-Gaussian statistics is summarized in the existence of 
anomalous scaling laws of velocity moments: 

{\UnY')-k-^^^\ (15) 

with C,{p) 7^ p/2C,{2). Anomalous scaling, also known as intermittency, is the quantitative way to state that velocity 
PDF's at different scales cannot be rescaled by any changing of variables. 



Let us now discuss two subtle points. Using some general arguments from the dynamical systems theory, one has 
that all the (typical) correlation functions at large time delay have to relax to zero with the same characteristic 
time, related to spectral properties of the Perron-Frobenius operator. If one uses this argument in a blind way, 
the apparently paradoxical result is that all correlation functions, C„^„(t) — (w„(t)M„(0)), must go to zero with same 
characteristic times. On the contrary one expects a whole hierarchy of characteristic times distinguishing the behavior 
of the correlation functions at different scales [|7|. In particular, the self-correlation fmiction, C„^„(i), decays with a 
characteristic time decreasing with n. The paradox is only apparent since the dynamical systems argument is valid at 
very long times, i.e. much longer than the longest characteristic time, and therefore in systems with many different 
time-fluctuations it is not helpful. In fact, it is well established numerically and well understood theoretically p7|-p^] 
that general multi-scale multi-time correlation functions of the kind C^'m(0 = (l'^n(0)l''l'"rn,(i)|'^) are described by the 
cascade formalism. In particular, most of the statistical properties in the inertial range can be well parametrized by 
the multi-fractal-formalism. 

On the other hand, the response properties are related to infinitesimal perturbations. Therefore, a priori, it is not 
obvious that the response depend on inertial range properties. The existence of the F/R relation and the fact that 
Cn.n{t) (and other similar correlation functions) are determined by the inertial range properties suggest that also the 
response features are ruled by the inertial range behavior if the invariant measure is dominated by local interactions 
among shells. 

Let us now examine the numerical results concerning Response functions in the shell model. 

Figure |j shows the diagonal mean response, GJJ(i), for a range of inertial shells, n S [7—14]. The most striking 
property is the impossibility to follow the response behavior at large scales (small shells) for large times, i.e. the 
explicit evidence that errors grow exponentially. In order to compare the different behavior (and different error 
propagation) between response and correlation function, we plot in Fig. (ra) both the average response and the self- 
correlation, Cn,n{t), for the shell n = 10. As it is clear from the previous figures, only response at the smallest scales 
(fast scales) in the inertial range can be computed with enough accuracy to follow an asymptotic decay. Still, also for 
this response the clear departure between the response and the self-correlation shows is another indication that the 
inertial-range statistics is far from Gaussian. By using HTS we can get an information on the temporal dependence of 
response function in the whole range of inertial shells. In Figure ^, we report some realizations of the instantaneous 
response function, which shows typical halving time experiments. 

In Figure (R), we summarize the results we obtain by comparing the mean halving time, {T{n)), with the character- 
istic times one extract from the decay properties of both mean response, Tij(n), and correlation functions, Tc{n), for 
those shells where such a behavior can be safely extracted. It is worth noticing how the mean halving time allows a full 
characterization of time properties also for those shells where the mean response GJ^ (i) cannot be measured for large 
time legs, t. Also, the dependence from the scale of the mean halving time is given as a best fit (T(n)) ^ k~^, with 
X = 0.53 ± 0.03. The value x = 0.53 ± 0.03 can be seen as an intermittent correction to the dimensional inertial-range 
prediction 2/3. On the other hand, the dependency from the scale of Tji{n) is difficult to extract due to the small 
number of points available. 

Let us now focus on the whole PDF of the halving time statistics. We first analyze the positive and negative 
moments of the halving times: 

TP{n) ^ {{TP{n)) ^ k-^^P\ (16) 

with p = — 5, ..., 3. Dimensional, non intermittent, scaling would predict the linear behavior for the scaling exponents: 
4'{p) = 2/3p. In Figure p^ we plot the results for the halving times scaling exponents 4'{p) ^or all moments from p = 
—5, . . . , 3 and the straight line corresponding to the dimensional inertial range estimate. We notice that intermittent 
corrections are much stronger for the positive moments than for the negative moments. This must be related to the 
fact that positive moments of the halving times are dominated by rare events where the response has a very long 
decaying. We interpret the fact that tlj{p) ~ (2/3)p as an indication that linear Response Functions are inertial range 
quantities. The latter results leads to the important conclusion that the invariant measure is well approximated by 
short range interaction among shells in the inertial range. 

As for the non-diagonal response function and for the generalized response function of higher order the numerical 
problems to measure them are even more pronounced. First, let us examine the off-diagonal response function 
R™{t) — Sun{t) / SumiO) . Of course these responses start from zero at time zero instead than from one as in the 
diagonal case. Measure them by a direct average is strictly impossible because of very large errors. We still decided 
to measure their characteristic time by using the time the response reach a "macroscopic" fraction, say 1/2, of the 
typical fluctuations on the scale where we are measuring the response, i.e. we collect the statistics of the first times r 
such that -R™(^ = r) = l/2(|u„p). We expect a strong asymmetry of the characteristic times depending whether the 
perturbation is done at smaller (m > n) or larger {n > m) scales. Indeed, by using usual inertial range arguments 
we expect that the response reacts always with a typical time given by the time of the largest between the two shells 



involved n, m. By fixing therefore the shell where we perturb, say m, we expect that the typical time of i?™(i), r„ ' 
is constant when n > m and scales as fc„ when n < m. In Figure ^ we show that indeed this behavior is well 
reproduced numerically. 

The strong intermittency shown by halving times in Figure ^ is the clear signature of deviations from simple Gaussian- 
like behavior of response functions. We must therefore also expect that generalized responses of higher order are not 
simply related to the linear response. For example, let us consider the third moments of the linear response [P0| : 

A simple non-intermittent behavior would suggest that Sl^{t) oc {G'^{t))^ while in Figure [lO| it is possible to see that 
this is definitely not the case for all t legs where we have a measurable signal. Unfortunately the already discussed 
statistical problems in measuring averaged response functions for long times are even more pronounced for generalized 
response functions. Therefore we refrain from showing any results for the scaling behavior of typical times of the 
generalized response functions. 



III. CONCLUSIONS 

We have addressed the problem of measuring time-properties of response functions in Gaussian models (Orszag- 
McLaughin) and strongly non-Gaussian models like shell models for turbulence. We have introduced the concept 
of halving time statistics with the aim to have a statistically stable tool to quantify the time decaying of response 
functions and generalized response functions of high-order. We have shown numerically that in shell models for three 
dimensional turbulence Response functions 

are inertial range quantities. This is a strong indication that the invariant measure describing the shell-velocity 
fluctuations is characterized by short range interactions between neighboring shells. 

Response functions and generalized response-functions play an important role in any diagrammatic approach for non- 
linear out-of-equilibrium systems. In this work we have presented the first numerical attempt to measure some of their 
properties in a systematic way. More work is needed, both numerical and analytical, in order to better understand 
the detailed structure of the invariant measure governing F/R relations. 

We acknowledge useful discussions with G. Boffetta and V. L'Vov. This work has been partially supported by the 
EU under the Grant No. HPRN-CT 2000-00162 "Non Ideal Turbulence". GL thanks the Department of Physics of 
L'Aquila University for the kind hospitality. 
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FIG. 1. Plot of the averaged Response Functions GZ(t) and Correlation Functions Cn,n{t) for five fast variables of the 
modified Orszag-McLaughlin model, n = 6 (+), n = 7 (x), n — 8 (*), n = 9 (□) and n — 10 (o). Statistical error bars are 
shown only for Response Functions corresponding to n = 6 and n = 10. Thin lines represent Correlation Functions. The 
statistics is over 10'' events. 
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FIG. 2. Log-log plot of Correlation Times Tc(n) (A), Mean Halving Times {r^n)) (D), and Response Times tr^u) (o) 
as function of k„, for the modified Orszag-McLaughlin model. Notice the much larger errors found when measuring the 
characteristic Response Times, Tnin). Errors on Tciji) and {r{n)) are of the same size as the representative symbols. The 
statistics is over 10 realizations. All these characteristic times follow the same scaling law with fc„. The exponent —3/2 of the 
scaling law follows from dimensional arguments. 
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FIG. 3. Collapse of the rescaled PDFs of the halving times for the modified Orszag-McLaughlin model. For simplicity, 
only the PDFs relative to the fastest four variables are shown, n = 7, ..., 10. The statistics is over lO"" impulsive infinitesimal 
perturbations as in Figure 0. 
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FIG. 4. Modulus of the average response functions, G^{t) — {R^{t)), for shells n — 7, ... ,14 (from top to bottom). Error 
bars are shown only for the smallest and the largest scales. The number of independent kicks used to perform the averages is 
around 2.10^. Notice the extremely large error bars measured for the slowest shell variables. The parameters entering in the 
eqs. of motion (Q are b = 0.4, f = 5.10"^, for iV = 25 shells. 




t 
FIG. 5. Comparison between the averaged Response Function, GJl(i), (top) and the self-correlation, Cn.n{t) (bottom) for 
the shell n = 10. Notice the different order of magnitude of error bars. 
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FIG. 6. Plot of three different diagonal instantaneous responses, R]^{t), for the shell n = 10, versus time. Halving time is 
fixed by the first time when the curve touches the threshold at A = 1/2. 
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FIG. 7. Log-log plot of the mean halving times, {T(n)) (+) and of decaying time of the mean diagonal response, Tii{n) (x), 
versus k„. We have checked that a different choice of the threshold A — 1/2 used to compute halving time does not affect the 
slope of the graph. 
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FIG. 8. Vp exponents of the p-th moment of halving times. The straight Une corresponds to the dimensional prediction 
V'p = 2/3p. 
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FIG. 9. Log- log plot of the off-diagonal response characteristic times r„ ' versus A:„. We performed two experiments. First 
we perturb at large scales, m = 2, and we follow the response at small scales n > 2, (-I-); the expected independence of 
characteristic times from the scale is well reproduced. Second, we perturb at small scales, m = 13, and we follow the response 
at larger scales, n < 13, x . In the latter case, for comparison we also plot the straight line (dashed) with the expected 
dimensional slope —2/3. 



13 




FIG. 10. Comparison between the third power of the mean diagonal response, 
response, \Sn{t)\, (x) computed for the shell n = 13, with 5.10^ kicks. 



IGJKi)!"^, (+) and the generahzed third order 
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